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The radio approach for detecting the ultra-high energy cosmic neutrinos has become a mature 
field. The Cherenkov signals in radio detection are originated from the charge excess of particle 
showers due to Askaryan effect. The conventional way of calculating the Cherenkov pulses by 
making Fraunhofer approximation fails when the sizes of the elongated showers become comparable 
with the detection distances. We present a calculation method of Cherenkov pulses based on the 
finite-difference time-domain (FDTD) method, and attain a satisfying effeciency via the GPU- 
acceleration. Our method provides a straightforward way of the near field calculation, which would 
be important for ultra high energy particle showers, especailly the electromagnetic showers induced 
by the high energy leptons produced in the neutrino charge current interactions. 

PACS numbers: 



I. INTRODUCTION 

Cosmic neutrinos, as a probe of the universe to the 
highest energy regime, are wonderful in many respects. 
Due to their extremely small interaction cross section, 
they can penetrate through galactic infrared (IR) and 
cosmic microwave background (CMB) photons, while 
photons of energy above 10 TeV would be attenuated. 
Furthermore, being uncharged, they propagate along 
straight lines and therefore are able to point directly back 
to their sources, while protons or other charged particles 
would be deflected by the magnetic field in the universe. 

Ultra-high energy cosmic rays (UHECRs) have been 
observed up to ~ 10 19,6 eV. The source of such amazingly 
energetic events have remained a mystery. Above this en- 
ergy scale, UHECRs interact with CMB photons through 
the Greisen-Zatsepin-Kuzmin(GZK) processes pQ. The 
GZK cut-off of the cosmic ray energy spectrum has been 
first observed by the High Resolution Fly's Eye Exper- 
iment [2] an d later confirmed by the Pierre Auger Ob- 
servatory [3], so the corresponding GZK neutrinos are 
almost guaranteed to exist. Nevertheless, none of these 
have been observed so far. Detecting the GZK neutrinos 
provides critical informations for unraveling the mystery 
of the origin and evolution of the cosmic accelerators, and 
will be one of the utmost tasks in the coming decade [4] . 

One promising way of detecting UHE neutrinos is the 
radio approach. When an ultra-high energy cosmic neu- 
trino interacts with ordinary matters on the Earth, it 
would lead to a hadronic debris, either by charged cur- 
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rent or neutral current. The former also produces a lep- 
ton with corresponding flavor. Both the high energy lep- 
tons and the hadronic debris induce particle showers. As 
proposed by Askaryan in the 1960's [5], the high energy 
particle shower develops in a dense medium would have 
net negative charges. This charge imbalance appears as 
a result of the knocked-off electrons being part of the 
shower, as well as the positrons in the shower annihilat- 
ing with the electrons of the medium. The net charges of 
the showers, typically 20% of total shower particles, serve 
as a source emitting the Cherenkov radiations when they 
travel in the medium. The sizes of the showers are quite 
localized (tens of cm in radial and few meters in longi- 
tudinal development) compared to those develope in the 
air (km scale), and therefore result in coherent radiations 
for wavelengths longer than the shower sizes. The cor- 
responding coherent wavelength turns out to be in the 
radio band, from hundreds of MHz to few GHz. This 
Askaryan effect has been confirmed in a series of exper- 
iments at Stanford Linear Accelerator Center (SLAC), 
using different dense media such as silica sand, rock salt 
and ice jBH9]. 



II. COHERENT CHERENKOV PULSES 

In the radio detection experiment, the signals come 
from the Cherenkov radiations of the net charges in the 
shower. The key concept which makes this technique 
possible is the coherent emission. In fact, the Cherenkov 
radiation is a broad band emission and the intensity in- 
creases as frequency. For a single charged particle, the 
Cherenkov signal in the radio band should be the weak- 
est in the spectrum. It is the compact size of the shower 
that makes radio signal so special. The coherent emission 
greatly enhences the signal strength in the radio band. 
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The electric field of Cherenkov radiations can be calcu- 
lated by solving the inhomogeneous Maxwell equations, 
as it has been demonstrated in the paper of Zas, Halzen, 
and Stanev [10] . The vector potential can be obtained 
by the Green's function method: 

exp (ik\x — x'\) 
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where k is the wavenumber, x the position of the de- 
tector, x ' the position of the shower particles, and J the 
current sources. Adopting the Fraunhoffer approxima- 
tion 

exp (ik\x — x ; |) exp (ik\x\ — ikx • x ') , , 

|x-x'| ~ S ' ( ] 

where R is the absolute value of x, the integration in 
Eq. ( [T]) can be greatly simplified and therefore enhances 
the computational effeciency in the Monte Carlo simula- 
tion [TQHT21 [T7] . The validity of this approximation relies 
on several length scales: the detection distance (R), the 
spatial size of the shower (I) and the wavelengths of in- 
terest (A). The Fraunhoffer approximation works well 
under the condition 
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where is the observational angle between the shower 
axis and the observational direction. 

However, for the ultra-high energy showers the lon- 
gitudinal development is longer, especially for the 
electromagnetic showers that suffer from the Landau- 
Pomeranchuk-Migdal (LPM) suppression [T3l [14] . Elec- 
tromagnetic showers can be produced by the charge cur- 
rent generated leptons. For a electromagnetic shower 
of EeV-scale energy, The impact of LPM effect on the 
shower development has been investigated by Monte 
Carlo simulations [T5HT7] . Electromagnetic showers of 
primary energies 10 20 eV can be extended to about 200- 
m long with great fluctuations. In such cases, the far 
field condition cannot be satisfied for distance up to sev- 
eral kilometers, while the typical detection distance for 
ground array detectors is about 1 km due to the attenu- 
ation length of radio signals in ice. Under these circum- 
stances, the Fraunhoffer approximation is clearly invalid, 
and one has to deal with the complicated integration in 
Eq.(0. 

In the paper of Buniy and Ralston [18 , the correction 
has been made by the saddle point approximation, while 
it still cannot cope with the extreme cases for R ~ Z, the 
near field regime. We handle this problem by a numerical 
method based on first principle so that the near field 
radiations can effeciently obtained. Although far field 
radiations would be more time consuming, it is not our 
focus in this paper. 

Hadronic showers, on the other hand, are less affected 
by LPM effect JTSJ [20] , since the sources of electromag- 
netic components in hadronic showers are the decay of 
neutral pions, which tend to be interact with matters in- 
stead of decay at energy above 6.7 PeV. The far field 
condition is well fulfilled for hadronic showers in most 
parctical cases. 



III. NUMERICAL METHOD 

Numerical algorithms for calculating electromagnetic 
fields have been existed for decades. However, it was 
not until the recent rapid growth of computational power 
that this approach became wildly adopted. Among all 
the existing algorithms, the finite- difference time-domain 
(FDTD hereafter) method has several advantages: 

• It is exceptionally simple to be implemented by 
computer programs. 

• It is a time-domain approach that is well suitable 
for an impulse signal. A broad band impulse can 
be calculated in one single run. 

• The algorithm itself is inherently parallel and the 
effeciency can be largely improved via parellel com- 
puting. 

The idea of FDTD was first proposed by Yee in the 
1960's [21], and has been in use for many years for 
the electromagnetic impulse modeling. Like most of the 
numerical finite difference methods, the space are dis- 
cretized into small grids and fields are calculated on each 
grid by solving Maxwell equations. Adopting a special 
lattice arrangment (known as the Yee lattice), the E-field 
and H-field are staggered in both space and time and can 
be calculated in a leapfrog time-marching way. 



A. Algorithm 

The Maxwell curl equations in differential forms are 

dU 



V x E = -/i 



dt 



V x H = e^-+aE 
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The FDTD method approximates derivatives by finite 
differences. The central difference is adopted to achieve 
2nd order accuracy in both spatial and temporal deriva- 
tives. We assume cylindrical symmetry along the shower 
axis (defined as z-axis), and therefore all the derivatives 
with respect to <p vanish. In addition, due to the polar- 
ization property of Cherenkov radiations, the H ri H Z1 
components also vanish. This can save large amounts 
of computer memories as well as calculation time. Fig- 
ure ([!]) shows the configuration of the lattice under these 
assumptions. Maxwell equtions in a cylindrical coordi- 
nate are reduced as: 
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We can use Eq.( (J6|) and Eq.( Q) to update the E- 
field, and then use Eq.( (fSJ) ) to update the H-field. The 
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spatial and temporal grid sizes have been chosen such 
that the numerical stability is satisfied and the numerical 
dispersion is controlled at an acceptable level [22) [23]. 
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FIG. 1: The lattice configuration we adopted. 



B. Simulation Setup 



Cherenkov radiations. We define a spherical coordinate 
whose origin lies on the intersection of the z axis and the 
shower maximum. Choosing the desired detector posi- 
tions by varying detection angle (9) and detection dis- 
tance (R), and record the electric field values as time 
evolves. After one single run, we can have all the simula- 
tion results we want. Here one can see the merit of using 
a time-domain calculation method. 




In order not to lose the focus on the RF calculation, we 
simply assume a shower model for the longitudinal devel- 
opment rather than invoking the Monte Carlo packages. 
For electromagnetic showers, the well known Nishimura- 
Kamata-Greisen (NKG) parametrization formula [24] de- 
scribing the number of shower particles is 



N e (X) 
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where Eq is the energy of the primary particle and E c 
the critical energy, X the slant depth, X max depth as 
shower reaching its maximum number, s the (dimension- 
less) shower age defined as s = 3X/(X + X max ). The 
NKG formula can be fit by a Gaussian distribution as: 



N e (z) = iV max exp(- — ) 



(10) 



where iV max is the particle number at shower maximum 
and I is the longitudinal shower length. The NKG for- 
mula for Eq = 10 TeV has I ~ 2 m in ice. For charge 
distribution in a snapshot, we assume the Gaussian dis- 
tribution in both radial (r) and longitudinal (z) direc- 
tions: 



n(z, r) — 
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where z and r are the cylindrical coordinates in the unit 
of g/cm 2 , and a z and o r are the standard deviation of the 
distribution in z and r direction respectively. We choose 
<j r = 5 cm in our simulations according to the Moliere 
radius in ice. The a z is generally even smaller than the 
Moliere radius and thus has no significant effect on the 
radiation pattern. 

Figire ([2| is a cartoon that demonstrates the setup. 
The shower travels in the +z direction emitting 



FIG. 2: Cartoon for our simulation set up. 



IV. PERFORMANCE IMPROVEMENT VIA 
THE GPU PARALLELIZATION 

Recently, graphical processing units (GPU) have be- 
came more and more important in the field of high per- 
formance computation. Under the needs in the computer 
game markets, the GPUs are currently under booming 
developments. The GPUs have become programmable 
via the Compute Unified Device Architecture (CUD A) 
provided by NVIDIA, and have been implemented in var- 
ious scientific areas such as molecular dynamics, gravi- 
tational N-body simulations and lattice QCD [25H29] . 
CUD A is an extension of C language, one of the most 
popular high level languages in the world. Programmer 
familiar with C language can utilize GPU computations 
by simply calling the functions from CUD A. The gen- 
eral parallelization strategies can be found in the CUDA 
Programming Guide (available on their webpage). 

Because of the multicore architecture, GPUs are 
ideal for implementing parallel algorithms. The FDTD 
method is therefore the ideal candidate to benefit from 
GPU computing. The texture memories provide the pos- 
sibility to realize the memory cache speedup. Since the 
texture memories in CUDA are read-only, we bind the ID 
cuda array to E-field and H-field alternatively, and write 
back to global memories for the unbound E/H-field. The 
basic steps are like the following: 

1. Allocate global memories to store E-field and H- 
field. 
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2. Bind texture to H-field. 

3. Calculate and store the updated E- field by reading 
the cached H-field. 

4. Unbind H-field. 

5. Bind texture to E-field. 

6. Calculate and store the updated H-field by reading 
the cached E-field. 

7. Unbind E-field. 

8. Repeat steps (2) - (7). 

We have investigated the performances of our codes 
with one single NVIDIA GTX285 graphic card, which has 
240 cores and 933 floating point operations per second 
(GFIOPS) of theoretical peak performance. The follow- 
ing table shows the performances of our codes with differ- 
ent total grid numbers (N), compared to Intel Quad Core 
i7 920 at 2.66GHz (sequential code without CPU paral- 
lelization). The GPU system is provided by the Center 
for Quantum Science and Engineering of National Tai- 
wan University (CQSE). 



iVgrid 


tcpu 


£gpu 


speed up 


GPU performance 




(sec) 


(sec) 


(^cpuAgpu) 


(GFLOPS) 


256 2 


5.15 


0.040 


128.75 


57.0 


512 2 


34.94 


0.186 


187.85 


98.1 


1024 2 


294.41 


1.167 


252.28 


125.1 


2048 2 


2122.80 


8.372 


253.56 


139.5 


4096 2 


15751.17 


67.617 


232.95 


138.2 


8192 2 


120617.48 


681.082 


177.11 


109.8 



TABLE I: Performance of CPU and GPU codes. 

The algorithm of FDTD method is inherently paral- 
lel since each grid can be updated independently, so the 
advantage of GPU systems can be maximally utilized in 
this algorithm. The performance test shows a tremen- 
dous speed up using GPU parallelization, which allows 
the calculation to be done in very reasonable time, with 
a cost-efficient computer resource. 

V. RESULTS 

A. At the Cherenkov Angle 

First we compare the E-fields fixed at Cherenkov an- 
gle (9 C ) with various distances (R). The magnitude of 
the E u decreases as R increases. However, the decreas- 
ing speeds in different frequencies are not the same. For 
example, Fig. ([3]) shows the E u spectra at R = 25 m and 
R = 50 m, but with the latter multiplied by a factor of 
two. These two spectra match well at high frequencies 



while they deviate from each other at low frequencies. In 
Fig.( [4|, the spectra at R = 25 m and R = 50 m are 
shown with the latter scaled by a factor of two. In this 
case, the two spectra match well at low frequencies, while 
they deviate at high frequencies. According to these two 
examples, one can see that the behavior in high frequen- 
cies suggests that E u oc 1/ y/R, namely a cylindrical wave 
behavior. On the other hand, the behavior in low fre- 
quencies suggests E u oc 1/R, & spherical wave behavior. 



The different behaviors in different frequency regimes 
can be interpreted in a physical way. It is a known fact 
that waves of higher frequencies are less likely to diffract. 
Therefore, the high frequency waves in Cherenkov radi- 
ations are more confined in the ^-direction, and their 
energies can only spread into the Cherenkov cone. Their 
energies go like 1/R due to geometrical reason, and hence 
1 / y/R for the fields. For the low frequency waves, diffrac- 
tion allows another direction (the ^-direction) for their 
energies to disperse, and therefore decrease faster. 



The fact that the higher frequency regime decreases 
slower implies that there is a shift of the peak frequecy 
in different R. The peak will migrate to the higher 
frequency regime as distance moves further away. The 
shape of the spectrum is independent, and a simple scal- 
ing relation of R is no longer valid here. 




FIG. 3: Eu spectra at R = 25 m and R = 50 m with the 
latter scaled by a factor of y/2\ these two spectra match 
well at high frequencies, implying a cylindrical behavior 
(E u oc 1/VR). 
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FIG. 4: E u spectra at R = 25 m and R = 50 m with the 
latter scaled by a factor of 2; these two spectra match well 
at low frequencies, implying a spherical behavior (E u ex 

w 



FIG. 5: The E„ - R relation of 100MHz. The blue curve 
is the E-field. The green and red curve show the cylin- 
drical (Eu ex 1/ \fR) and spherical {E u oc 1/R) behavoir 
respectively. The E-field goes like 1 / \fR in the near field 
regime, while approaches to (1/i?) -behavoir in the far 
field regime. Note the smooth transition from near field 
to far field. 



B. Angular Distribution 



In principle, all waves will eventually diffract and be- 
have like spherical waves no matter how high their fre- 
quencies are, if we set the distance R infinitely large. 
Therefore, the terms "high" and "low" frequencies are 
only relative concepts. From the far field condition in 
Eq.( |3| we can see that all three length scales couple 
with each other. The waves start to diffract after they 
propagate to the distance large enough such that the far 
field condition is satisfied. This character can be demon- 
strated more clearly if we plot the E u - R relation with 
one single frequency. Figure ([5| shows how E u decreases 
with R. At the distance very close to the shower, E u goes 
like 1/y/R. As the distance increases, for large enough i?, 
the radiation can be viewed as a point source and thus 
have E u oc 1/R. We can see a smooth transition from 
cylindrical behavoir to spherical behavoir. At R large 
enough such that all the waves of different frequencies in 
the Cherenkov pulse reach the far field regime, the shape 
of the spectrum is fixed and independent of i?, and R 
becomes just a normalization factor. 



The diffraction effect can also be seen in the angular 
distribution. Figures ( |6| and ( [7]) show the angular dis- 
tributions of E u at R = 4 m and R = 64 m respectively. 
For the case at R = 4 m, the distance is too close for 
the waves to diffract. In fact, the radiation pattern in 
the near field (before diffraction happen) is just a fuzzy 
image of the radiating source. We can see the waves of 
higher frequencies have stronger magnitudes, which is the 
inherent character of Cherenkov radiations. Note the dis- 
tribution is not symmetric on two sides, since the 6 < 6 C 
part is closer to the shower axis than the > C part is. 
This asymmetry can be understood as we are using the 
spherical coordinate to describe a sysytem that is actu- 
ally cylindrical symmetric. For the case at R = 64 m, the 
detection distance is long enough for waves diffract into 
the ^-direction. The lower the frequency is, the wilder 
of its angular distribution is, which is the standard prop- 
erty in diffraction. In time domain, it can be seen in 
Figure ( [8| that the pulse at detection angle more away 
from C has wilder width. In principle, if the frequency 
goes to infinity, the angular distribution will be a delta 
function. However the destructively interference of the 
lateral distribution suppress these high frequecy compo- 
nents. The distribution now looks much more symmetric, 
since the differences of the distances to the shower axis 
between 6 < 6 C and 6 > 6 C parts are negligible compares 
to R. Namely, as R goes further, the system, originally 
cylindrical symmetric, becomes more and more spherical 
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symmetric. 
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FIG. 6: E u angular distributions for R = 4 m (near field). 
The diffraction has not yet happened and the radiation 
pattern follows the shower image. 
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FIG. 8: Time domain E- fields at different angles at R 
= 32 m. The pulse is wilder at angles more away from 
# c , implying the fact that lower frequency components 
diffract more. 



C. Comparison 
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FIG. 7: E u angular distributions for R = 64 m (far field). 
The higher the frequency is, the narrower of its angular 
distribution is. Diffraction phenomenon is quite obvious 
here. 



We compare our results with the conventional far field 
formula. The one dimensional approach in [11] should 
be a reasonable approximation except at the Cherenkov 
angle. Substituting our shower model in Eq.( [To|, the 
E-field can be obtained: 
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where the parameter p(Q,uj) = (1 — ncosQ) uj/c. Figure 
( |9| shows the comparison between the spectra of the 
far field formula and our simulation results at R = 64 m 
and different 6. At lower frequency part they are in good 
agreements, while at high frequency part there are sig- 
nificant differences between them. It can be understood 
as the disagreement part has not yet reached the far field 
regime and thus decreases slower. If we set larger R, the 
disagreement part will enter the far field regime and will 
match with the formula. 
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FIG. 9: Comparison between the spectra of the far field 
formula and our simulation results at R = 64 m. From 
top to bottom are the spaetra at # c + 5°, C + 1O°, # c + 15°, 
C + 20° respectively. The solid curves are our results and 
the dashed curves are obtained from the far field formula. 



VI. SUMMARY AND CONCLUSIONS 

We have developed a numerical code to calculate the 
radiation patterns of Cherenkov signals from near field 
to far field based on the FDTD method. By utilizing 
GPU parallel computation, the effeciency of the code can 
be greatly improved to a satisfying level on a comercial 
graphic card NVIDIA GTX285. This will be useful in 
studying the signals originate from an elongated shower 
of its size comparable to the detection distance, where 
the traditional Fraunhofer approximation does not ap- 
ply. Signals from the ultra-high energy electromagnetic 
showers induced by the electrons produced in neutrino 
charged-current interaction are the typical examples, for 
they suffer from severe impact of the LPM effect. Our 
result shows a smooth transition between near field and 
far field pattern. In fact, the FDTD method is more 
suitable for the calculation of near field pattern since the 
far field pattern may challange the computer resources. 
The spectrum and the angular distribution of near field 
pattern have quite complicated R dependences instead of 
simple scalings in the case of far field. 

In the cases of far field, the angular distribution of sig- 
nals induced by LPM-elongated showers are much nar- 



rower than the ordinary ones, and the detection solid 
angle is considered to be small. However the far field as- 
sumption neglects the shower size and treat it as a point 
source which is not fair. A shower of hundred meters 
long would in fact generates signals spaneed also hun- 
dred meters, and is surely as possible to be detected as 
those with compact size. 

The idea of using staggered grid configuration to solve 
the two coupled first-ordered PDEs is not limited to 
the electromagnetic problems. Recently there are also 
applications of FDTD method in the acoustic simula- 
tions [30j El] solving the coupled pressure fields and ve- 
locity fields. It is possible to simulate signals in the neu- 
trino detection experiments using acoustic approaches, 
which is another potential field in UHE neutrino detec- 
tion [32j[33]. 

In the next generation neutrino detectors applying the 
ground array layout, it is possible to simultaneously de- 
tect the hadronic shower and the electromagnetic shower 
that are induced by one single charged current neutrino 
interaction. If both the hadronic shower and the electro- 
magnetic shower can be correctly reconstructed, it opens 
an opportunity to distinguish the electron neutrino from 
others [34], since the two shower vertexes are nearly lo- 
cated at the same places. However, the features of near 
field are very different from far field, and will face some 
detection difficulties. For example, the normal way to 
reconstruct the direction of incoming Cherenkov pulses 
by the arrival time differences between antennas is based 
on the assumption that the shower is a point source, i.e. 
the far field assumption. For a extended shower this as- 
sumption fails and therefore requires a new reconstruc- 
tion method. Any ground based neutrino detector has to 
take this near field effect into account in order to recon- 
struct signals from extended showers. 
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